% To do likelihood ratio test
clear all;
clc;
disp('To do likelihood ratio test_MK3');

%load MK2 data (Both tests based on same data)
% load PathCountTotal_MK3_TP06YRJfinR5G1000kMR7.mat;
% load PathCountAvgProb_MK3_TP06YRJfinR5G1000kMR7.mat;
load ('MK3_TP06YRJfinR5G2501000kMR7.mat','PathCountTotal','PathCountAvgProb');
filename='TP06_Y_RJ_fin_R5_G250_1000k_MR7_Likelihood Ratio Test_MK2_MK3.xls';

% take the avg data from PathCountAvgProb (MK2).
N3(:,:,:,:) = PathCountTotal(1,:,:,:,:);
P3(:,:,:,:) = PathCountAvgProb(1,:,:,:,:);
ttlData = sum(N3(:));

% to calculate log (LMK2) hyphothesis
N2(:,:,:) = sum(N3(1:7,:,:,:));
P2=zeros(7,7,8);
LLMK2=0;
for i=1:7
    for j=1:7
        for k=1:8
        P2(i,j,k)=N2(i,j,k)./sum(N2(i,j,1:8));
        if N2(i,j,k)~=0
            Y =log(P2(i,j,k)).* N2(i,j,k);
            LLMK2=LLMK2+Y;
        else
            LLMK2=LLMK2+0;
        end
        end
    end
end


% to calculate log (LMK3) hyphothesis
LLMK3=0;
for i=1:7
    for j=1:7
        for k=1:7
            for d=1:8                
                if N3(i,j,k,d)~=0
                    YY =log(P3(i,j,k,d)).* N3(i,j,k,d);
                    LLMK3=LLMK3+YY;
                else
                    LLMK3=LLMK3+0;
                end
            end
        end
    end
end


%Likelihood Ratio Test
invalid_Num_N2 = prod(size(find(N2==0)));
invalid_Num_N3 = prod(size(find(N3==0)));
Num_N2=prod(size(P2));
Num_N3=prod(size(P3));
dof = (Num_N3 - invalid_Num_N3) - (Num_N2 - invalid_Num_N2);
% dof = Num_N3 - Num_N2 - invalid_Num_N3 - invalid_Num_N2;% minus the number of zero element in P2 and P3.
[LRh,LRp,LRstat,cV] = lratiotest(LLMK3,LLMK2,dof);


% PP2(:,:)=P2(2,:,:);
% NN2(:,:)=N2(2,:,:); 
% PP3(:,:)=P3(2,2,:,:);
% NN3(:,:)=N3(2,2,:,:); 
% sheetlist1=strcat('Prob of ijk_pathR2');
% sheetlist2=strcat('Prob of ijkd_pathR2R2');
% sheetlist3=strcat('ObservationNum N2_pathR2');
% sheetlist4=strcat('ObservationNum N3_pathR2R2');
sheetlist5=strcat('Parameters');
% xlswrite(filename, PP2, sheetlist1);
% xlswrite(filename, PP3, sheetlist2);
% xlswrite(filename, NN2, sheetlist3);
% xlswrite(filename, NN3, sheetlist4);
d = {'LRh','LRp','LRstat','cV','dof','LLMK2','LLMK3','ttlData'; LRh, LRp, LRstat,cV,dof,LLMK2,LLMK3,ttlData};
xlswrite(filename, d, sheetlist5)


